Portfolio basics

This is the companion notebook to the introductory course on simple portfolio construction and performance metrics. It contains simple examples of portfolio backtesting. The structures of this section will be used later on in the course.

\(\rightarrow\) ABOUT THE NOTEBOOK:

• the code chunks are sequential. They must be executed in the correct order.
  
• text between stars appears in bold in html format.

1 A glimpse at the S&P500

The S&P500 is one of the most widely scrutinised index in the US Equity investment space. It serves very often as benchmark. The purpose of this section is to unveil a few of its statistical properties and to use the highcharter package for highchart rendering in R.

1.1 Prices

We start by loading the packages and downloading the data. For technical reasons, we download the series of the ETF replicating the S&P500 (SPY ticker).

if(!require(highcharter)){install.packages("highcharter")} # Package for nice financial graphs
if(!require(scales)){install.packages("scales")}           # Package for graph scales
library(tidyverse)   # The data wrangling package
library(plotly)      # For interactive graphs
library(lubridate)   # For data management
library(quantmod)    # The package that eases the downloading of financial data
library(highcharter) # The package for financial time-series plots

min_date <- "1980-01-01"
max_date <- "2021-05-05"
prices <- getSymbols("SPY", src = 'yahoo',  # Yahoo source (ETF ticker) ^GSPC
                     from = min_date, 
                     to = max_date,
                     auto.assign = TRUE, 
                     warnings = FALSE) %>%
  map(~Ad(get(.))) %>% 
  reduce(merge) %>%
  `colnames<-`("SPY")

Then, we turn to plotting, using highchart format. We underline that this package works with special xts (R extensible time-series) format. We point to the package reference for more details on this subject:
- https://jkunst.com/highcharter/
- https://www.highcharts.com/blog/tutorials/highcharts-for-r-users/
- https://cran.r-project.org/web/packages/highcharter/index.html

highchart(type = "stock") %>%
  hc_title(text = "Evolution of the SPY") %>%
  hc_add_series(prices)

A nice feature of highcharts is that they allow the user to change the observation period and see the values of points on the curve.

1.2 Returns

Next, we turn to the distribution of returns.

returns <- prices %>%                          # Formula for returns
  data.frame(Date = index(.)) %>%              # Adding the Date as a column
  mutate(SPY = SPY / dplyr::lag(SPY) - 1) %>%  # Returns
  na.omit()                                    # Removing NA rows
m <- mean(returns$SPY)                         # Average daily return
s <- sd(returns$SPY)                           # Volatility = sd of daily returns

returns %>% ggplot() +                         # Plot
  geom_histogram(aes(x = SPY, y = ..density..), bins = 100, alpha = 0.6) + theme_light() +
  stat_function(fun = dnorm, args = list(mean = m, sd = s), aes(color = "Gaussian")) +
  theme(legend.position = c(0.7, 0.7))

We plot the Gaussian density with parameters corresponding to the sample mean and standard deviation. Small grey rectangles around \(\pm 0.05\) indicate that large positive and negative returns occur more often than estimate by the Gaussian law: the tails of their distribution are notoriously heavy.

1.3 Volatility

Finally, we take a dynamic look at the volatility. We take a frugal approach; much more elegant methods are presented in Section 4.5 of Reproducible Finance by Jonathan Regenstein, as well as the roll_sd() function used in http://www.reproduciblefinance.com/drafts/post/vix-and-realized-volatility-updating-our-previous-work/

nb_days <- 63                                       # 63 days roughly equivalent to 3 months 
vol <- 0                                            # Initialisation
for(i in 1:(nrow(returns) - nb_days + 1)){          # Loop on dates: not elegant!
  vol[i] <- sd(returns$SPY[i:(i + nb_days - 1)])    # Forward vol computed on rolling window of nb_days
}
Date <- returns$Date[nb_days:nrow(returns)]         # Vector of dates
vol <- data.frame(vol*sqrt(252))                    # Vector inside a dataframe

rownames(vol) <- Date                               # Change the index of the dataframe
highchart(type = "stock") %>%                       # Code for the plot
  hc_title(text = "Evolution of realized volatility") %>%
  hc_add_series(as.xts(vol))

Clearly the graph shows periods of low volatility and clusters of high market turbulence (crashes, most of the time). Returns are therefore not stationary. The properties we exhibited for the S&P500 are also true at the individual stock level.

In the next sections, we turn to the core topic of portfolio backtesting.

2 Data

First, let’s start with the preprocessing of the data.

We first import and arrange the data. The data consists of monhtly financial information pertaining to 30 large US firms. They are characterised by their ticker symbol:

A - F	G - M	O - Z
AAPL (Apple)	GE (General Electric)	ORCL (Oracle)
BA (Boeing)	HD (Home Depot)	PFE (Pfizer)
BAC (Bank of America)	IBM	PG (Procter & Gamble)
C (Citigroup)	INTC (Intel)	T (AT&T)
CSCO (Cisco)	JNJ (Johnson & Johnson)	UNH (United Health)
CVS (CVS Health)	JPM (JP Morgan)	UPS
CVX (Chevron)	K (Kellogg)	VZ (Verizon)
D (Dominion Energy)	MCK (McKesson)	WFC (Wells Fargo)
DIS (Disney)	MRK (Merck)	WMT (Walmart)
F (Ford)	MSFT (Microsoft)	XOM (Exxon)

There are 7 attributes: closing price (Close), market capitalisation in M$ (Mkt_Cap), price-to-book ratio (P2B), 1 month volatility (Vol_1M), 1 month relative strength index (RSI_1M), debt-to-equity ratio (D2E) and profitability margin (Prof_Marg).
Finally, the time range is 2000-2021.

load("data.RData")                      # Loading the data: IF DIRECTORY OK!
data <- data %>% arrange(Date,Tick)     # Ranked first according to date and then stocks
summary(data)                           # Descriptive statistics

      Tick           Date                Close             Vol_1M           Mkt_Cap             P2B                 D2E          
 AAPL   : 255   Min.   :1999-12-31   Min.   :  0.217   Min.   :  5.867   Min.   :   4467   Min.   :   0.1013   Min.   :    0.00  
 BA     : 255   1st Qu.:2005-03-31   1st Qu.: 19.166   1st Qu.: 15.970   1st Qu.:  67392   1st Qu.:   1.4223   1st Qu.:   31.83  
 BAC    : 255   Median :2010-07-30   Median : 32.307   Median : 21.835   Median : 137516   Median :   2.3494   Median :   66.76  
 C      : 255   Mean   :2010-07-30   Mean   : 49.261   Mean   : 27.054   Mean   : 160615   Mean   :  11.1931   Mean   :  218.70  
 CSCO   : 255   3rd Qu.:2015-11-30   3rd Qu.: 57.307   3rd Qu.: 31.894   3rd Qu.: 210966   3rd Qu.:   4.2259   3rd Qu.:  213.09  
 CVS    : 255   Max.   :2021-02-26   Max.   :442.987   Max.   :268.563   Max.   :2255969   Max.   :1678.9053   Max.   :14585.61  
 (Other):6120                                                                                                                    
   Prof_Marg           ESG_rank      
 Min.   :-106.260   Min.   :  7.317  
 1st Qu.:   5.008   1st Qu.: 54.237  
 Median :  10.312   Median : 72.222  
 Mean   :  11.712   Mean   : 68.352  
 3rd Qu.:  19.171   3rd Qu.: 86.364  
 Max.   : 145.583   Max.   :100.000  
                    NA's   :5214     

tick <- unique(data$Tick)               # Set of assets
t_all <- unique(data$Date)              # Set of dates
N <- n_distinct(data$Tick)              # Number of stocks

This simple table shows a possible outlier for the P2B variable. The maximum value is clearly out of range.

Next, we format the data for future use. Notably, we compute & store returns.

data <- data %>% 
  group_by(Tick) %>%                                   # Grouping: returns computed stock-by-stock
  mutate(Return = Close / dplyr::lag(Close) - 1) %>%   # Adding returns
  ungroup()

returns <- data %>%                                    # Take data
  select(Tick, Date, Return) %>%                       # Select 3 columns
  pivot_wider(names_from = Tick, values_from = Return) # Put them into 'matrix' format
returns                                                # Show the returns

3 What students often do…

Now, let’s consider the usual Markowitz allocation subject to the budget constraint \(\textbf{w}'\textbf{1}=1\):

\[\text{max} \ \textbf{w}'\boldsymbol{\mu} - \frac{\gamma}{2}\textbf{w}'\boldsymbol{\Sigma}\textbf{w}, \quad \text{s.t.} \quad \textbf{w}'\textbf{1}=1 \] with Lagrangian: \[L(\textbf{w})=\textbf{w}'\boldsymbol{\mu} - \frac{\gamma}{2}\textbf{w}'\boldsymbol{\Sigma}\textbf{w}+\delta (\textbf{w}'\textbf{1}-1)\] so that \[\frac{\partial L}{\partial \textbf{w}}=\boldsymbol{\mu}-\gamma\boldsymbol{\Sigma}\textbf{w}+\delta\textbf{1}=0 \quad \Rightarrow \quad \textbf{w}=\gamma^{-1}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}+\delta \textbf{1}),\] where \(\delta\) is chosen to satisfy the budget constraint. A particular case is the minimum variance portfolio which corresponds to \(\boldsymbol{\mu}=\textbf{1}\): the optimal portfolio is proportional to \(\boldsymbol{\Sigma}^{-1} \textbf{1}\).

Sigma <- returns %>%                                       # Covariance matrix
  select(-Date) %>% 
  cov(use = "complete.obs")
Sigma[1:9,1:9] %>% round(3)                                # A quick look at Sigma

      AAPL    BA   BAC     C  CSCO   CVS   CVX     D   DIS
AAPL 0.013 0.002 0.003 0.004 0.004 0.001 0.002 0.000 0.003
BA   0.002 0.009 0.004 0.006 0.003 0.002 0.003 0.001 0.003
BAC  0.003 0.004 0.013 0.011 0.003 0.002 0.002 0.001 0.004
C    0.004 0.006 0.011 0.015 0.005 0.002 0.003 0.001 0.005
CSCO 0.004 0.003 0.003 0.005 0.009 0.001 0.002 0.000 0.003
CVS  0.001 0.002 0.002 0.002 0.001 0.006 0.002 0.001 0.002
CVX  0.002 0.003 0.002 0.003 0.002 0.002 0.004 0.001 0.002
D    0.000 0.001 0.001 0.001 0.000 0.001 0.001 0.002 0.001
DIS  0.003 0.003 0.004 0.005 0.003 0.002 0.002 0.001 0.006

Min_var_weights <- solve(Sigma) %*% rep(1,N)               # Proportional weights
Min_var_weights <- Min_var_weights / sum(Min_var_weights)  # Exact (scaled) weights
bind_cols(tick, Min_var_weights) %>%
  ggplot(aes(x = Min_var_weights, y = tick)) +             # Plot
  geom_col(alpha = 0.5) + theme_light()

Naturally, the stocks with super low volatility (Disney, Walmart) are favored! Let’s see when the weights are driven by the mean returns…

mu <- returns %>%                                    # Mean vector 
  select(-Date) %>%                                  # Remove date
  colMeans(na.rm = T)
Marko_weights <- solve(Sigma) %*% mu                 # Proportional weights
Marko_weights <- Marko_weights / sum(Marko_weights)  # Exact (scaled) weights
bind_cols(tick, Marko_weights) %>%
  ggplot(aes(x = Marko_weights, y = tick)) +         # Plot
  geom_col(alpha = 0.5) + theme_light()

AAPL and UNH have high returns => high weights.

Now let’s look at 2 ways of processing these weights into portfolio values.

w_matrix <- matrix(Marko_weights, 
                   ncol = N, nrow = length(t_all) - 1, byrow = T) # Weight matrix
(rowSums((returns %>% select(-Date)) * w_matrix) + 1) %>%         # Simple product = element-wise
  na.omit() %>%                                                   # Remove missing points
  cumprod() %>% round(3)                                          # Cumulative product + rounding

  [1]   0.957   0.895   0.931   0.933   0.927   0.996   0.985   1.097   1.085   1.095   1.101   1.209   1.236   1.256   1.300   1.359   1.310
 [18]   1.316   1.412   1.411   1.312   1.274   1.359   1.405   1.516   1.526   1.673   1.898   1.879   1.858   1.702   1.729   1.671   1.593
 [35]   1.538   1.607   1.639   1.564   1.597   1.630   1.781   1.846   1.871   1.886   1.835   1.928   1.859   1.954   1.960   2.043   2.138
 [52]   2.111   2.240   2.366   2.368   2.516   2.706   2.904   3.293   3.341   3.570   3.697   3.769   3.665   3.762   3.857   4.055   4.219
 [69]   4.486   4.718   4.957   5.056   5.044   5.060   4.884   4.965   4.759   4.843   5.204   5.171   5.360   5.560   5.720   5.700   5.829
 [86]   5.894   6.004   6.103   6.269   6.035   6.066   6.167   6.389   7.167   8.004   8.508   7.763   7.271   7.284   7.540   8.156   7.465
[103]   7.650   8.099   7.535   7.358   7.195   7.863   8.345   7.881   7.529   7.940   8.635   9.435   9.963   8.916   9.034   9.398  10.201
[120]  11.351  11.002  11.670  12.086  11.940  11.344  11.201  11.591  11.362  12.097  12.296  11.989  12.214  12.611  13.273  13.464  13.812
[137]  14.252  14.058  14.536  14.963  14.356  14.518  14.638  15.176  15.599  16.918  17.804  17.594  18.055  19.206  18.919  19.155  18.661
[154]  17.938  18.058  17.509  17.361  17.620  18.328  19.586  19.818  20.029  21.619  22.017  22.751  23.817  25.805  25.872  26.049  27.052
[171]  29.119  28.340  29.926  30.933  29.912  31.684  32.336  34.675  36.337  37.271  37.692  38.300  37.915  37.414  38.274  38.052  38.422
[188]  35.745  36.012  37.135  36.997  37.023  37.860  37.233  39.859  38.358  39.147  40.790  42.238  41.284  42.089  43.308  45.886  47.406
[205]  50.376  53.211  53.917  56.048  58.458  59.210  63.626  65.449  63.280  69.919  76.063  78.350  83.468  80.938  76.859  77.800  79.071
[222]  80.365  86.609  92.206  94.659  96.311 102.447  90.073  87.358  92.111  95.476  99.186  95.228  99.665  97.140 103.540 104.853 108.992
[239] 110.806 115.953 111.383  97.044  89.036 102.187 103.793  99.388  95.844 100.039  93.242  81.467  99.422  99.446  99.910  98.918

The above code is ugly. Because the data is neatly ordered, we can do better! Let’s add a new column to the data.

data %>% 
  select(Tick, Date, Return) %>%                          # Keep only necessary cols
  bind_cols(w_Marko = rep(Marko_weights, length(t_all)))  # Bind cols: beware to the order!

Now we can add benchmark weights and proceed.

p <- data %>% 
  select(Tick, Date, Return) %>%
  bind_cols(w_Marko = rep(Marko_weights, length(t_all))) %>%
  mutate(w_benchmark = 1/30) %>%                             # Add benchmark weights 1/N
  pivot_longer(cols = c(w_Marko, w_benchmark),               # Pivot to tidy
               names_to = "Portfolio",  
               values_to = "Weight") %>%
  group_by(Date, Portfolio) %>%                              # Group to compute portf return
  summarise(Return = sum(Weight * Return)) %>%               # Portf return, date-by-date
  ungroup() %>%                                              # Ungroup
  na.omit() %>%                                              # Remove outliers
  group_by(Portfolio) %>%                                    # Group along the 2 portfolios
  mutate(Value = cumprod(1+Return)) %>%                      # Compute portfolio values
  ggplot(aes(x = Date, y = Value, color = Portfolio)) +      # Plot
  geom_line() + theme_light() + scale_y_log10() +
  scale_color_manual(values = c("#706C66", "#0A5B90"), labels = c("Benchmark", "Markowitz")) + 
  theme(legend.position = c(0.3,0.7))

ggplotly(p, width = 900)

Given the log-scale: is this reasonably realistic? Obvious answer: NO! The weights have been chosen with knowledge of the trajectory!!! Also: portfolio managers update the weights…

4 A universe of functions

Backtesting portfolio strategies is usually simple at first, but becomes more intricate when many options are considered. To keep things simple, it is more efficient to work with functions. Functions help compartmentalise the different tasks of the process. We will build functions that compute portfolio weights and others that evaluate the performance metrics of the strategies. While we will usually work with a loop on backtesting dates, it is possible to build functions that directly generate portfolio returns (see below: the map() function).

By definition, a portfolio is a choice of weights that sum to one. Below, we implement two classical weighting schemes: the uniform portfolio (EW = equal weights) and the maximum Sharpe ratio portfolio (MSR). The latter is more complicated and requires inputs: namely the column vector of (expected) mean \(\mu\) and covariance matrix \(\Sigma\) of the assets. For simplicity, we will estimate them using sample moments - even though this is known to be a bad choice. This requires an additional argument in the function: the assets’ past returns. The MSR weights are \(w=\frac{\Sigma^{-1}\mu}{1'\Sigma^{-1}\mu}\). Both \(\mu\) and 1 are vectors here.

weights_msr <- function(returns){  # returns will refer to PAST returns
  m <- apply(returns, 2, mean)     # Vector of average returns
  sigma <- cov(returns)            # Covariance matrix
  w <- solve(sigma) %*% m          # Raw weights
  return(w / sum(w))               # Returns normalised weights
}

weights_ew <- function(returns){   # We keep the same syntax for simplicity
  N <- length(returns[1,])         # Number of assets
  return(rep(1/N,N))               # Weight = 1/N
}

We are now ready to proceed with the initialisation of the variables that will use in the main loop.

sep_date <- as.Date("2010-01-01")                # This date separates in-sample vs out-of-sample
t_oos <- t_all[t_all > sep_date]                 # Out-of-sample dates (i.e., testing set)
portf_weights <- matrix(0, nrow = length(t_oos), # Will store portfolio weights
                           ncol = length(tick)) 
portf_returns <- c()                             # Initialize portfolio returns

At last, we can proceed to the main loop. A note of caution: both in the weighting scheme functions and in the loop below, we assume well defined (finite, i.e., non NA) data. Obviously, feeding NA data in the system will produce NA outputs. There are only four steps in the loop:

1. extract the data
   
2. compute portfolio weights
   
3. compute realised returns
   
4. derive the return of the portfolio

for(t in 1:length(t_oos)){
  temp_data <- returns %>% filter(Date < t_oos[t]) # 1) Extracting past data: expanding window! 
  portf_weights[t,] <- temp_data %>%               # 2) Take this past data 
    select(-Date) %>%                                # Take out the date column
    na.omit() %>%                                    # Remove missing points
    weights_msr()                                    # Appply the weighting scheme function
  realised_returns <- returns %>%                  # 3) Take returns
    filter(Date ==  t_oos[t]) %>%                    # Keep only current date
    select(-Date)                                    # Take out the date column
  portf_returns[t] <- sum(portf_weights[t,] * realised_returns) # 4) Compute the return
}

In the above backtest, the amount of data considered to form the portfolio decision increases with time (expanding window). It is easy to fix the number of data point at each step and proceed on rolling windows (exercise). Likewise, switching from MSR to EW is immediate (just change the weighting function).

The main output is the vector of portfolio returns. We can easily plot the evolution of the portfolio through time.

port <- data.frame(Date = t_oos, 
                   Portfolio = cumprod(1+portf_returns))    # Portf. values via cumulative product
port %>% ggplot() + theme_light() +
  geom_line(aes(x = Date, y = Portfolio))                   # Plot

Below, we create a function that takes a series of returns as input and provides a few simple performance metrics.

perf_met <- function(returns){
  avg_ret <- mean(returns, na.rm = T)                     # Arithmetic mean 
  vol <- sd(returns, na.rm = T)                           # Volatility
  Sharpe_ratio <- avg_ret / vol                           # Sharpe ratio
  VaR_5 <- quantile(returns, 0.05)                        # Value-at-risk
  met <- data.frame(avg_ret, vol, Sharpe_ratio, VaR_5)    # Aggregation of all of this
  rownames(met) <- "metrics"
  return(met)
}
perf_met(portf_returns) # Let's test on the actual returns

Note: the values are not annualised. The annualisation can be directly coded into the function. A heuristic way to proceed is to multiply average returns by 12 and volatilities by \(\sqrt{12}\) - though this omits the compounding effect. These simplifications are ok if they are used to compare strategies. The Value at Risk is nonetheless dependent on the return horizon and cannot be proxied so simply.

An important indicator is the turnover of the portfolio: it assesses asset rotation and thus impacts transaction costs. Simple turnover computes average absolute variation in portfolio weights: \(\text{Turn}=\frac{1}{T}\sum_{t=2}^T\sum_{n=1}^N|w_t^n-w_{t-1}^n|\). Full turnover takes into account the variation of the weights between rebalancing dates: \(\text{Turn}=\frac{1}{T}\sum_{t=2}^T\sum_{n=1}^N|w_t^n-w_{t-}^n|\), where \(t-\) is the time just before the rebalancing date.

turnover_simple <- function(weights, t_oos){
  turn <- 0
  for(t in 2:length(t_oos)){                             # Loop on dates
    turn <- turn + sum(abs(weights[t,] - weights[t-1,])) # Variation in weights
  }
  return(turn/(length(t_oos)-1))                         # BEWARE: monthly value!
}
turnover_simple(portf_weights, t_oos)                    # Simple turnover!

[1] 0.2189094

Below, we switch to the full definition of the turnover. In this case, the variation in weights is adjusted because the weight prior to rebalancing have drifted.

turnover_full <- function(weights, asset_returns, t_oos){
  turn <- 0
  for(t in 2:length(t_oos)){
    realised_returns <- returns %>% filter(Date == t_oos[t]) %>% select(-Date)
    prior_weights <- weights[t-1,] * (1 + realised_returns)       # Before rebalancing
    turn <- turn + apply(abs(weights[t,] - prior_weights/sum(prior_weights)),1,sum)
  }
  return(turn/(length(t_oos)-1))
}
asset_returns <- returns %>% filter(Date > sep_date)    # Asset returns over the experiment
turnover_full(portf_weights, asset_returns, t_oos)      # Real turnover

[1] 0.3014361

Note that the turnover is computed at the monthly frequency. The second value is always more realistic; in this case it is substantially higher compared to the simplified proxy. For the sake of compactness, the turnover should be included in the perf_met() function.

5 Extensions

5.1 Comparing strategies

An important generalisation of the above framework is to consider more than one strategy. Comparisons are commonplace in the asset management industry: obviously, people look for the best strategy (according to particular goals, beliefs, and preferences). Below, we show how this can be handled. We will compare the two strategies that we mentionned above. First, we need to re-initiate the variables because their dimension will change (NOTE: an alternative route would be to work with lists).

Tt <- length(t_oos)                                             # Nb of computation dates 
# Avoid T because T = TRUE!
nb_port <- 2                                                    # Nb of portfolios
portf_weights <- array(0, dim = c(Tt, nb_port, length(tick)))   # Store weights
portf_returns <- matrix(0, nrow = Tt, ncol = nb_port)           # Store returns

Second, we embed all weighting schemes into one single function.

weights_multi <- function(returns, j){   # Strategies are indexed by j
  if(j == 1){ # j = 1 => MSR
    return(weights_msr(returns))
  }
  if(j == 2){ # j = 2 => EW
    N <- length(returns[1,])
    return(rep(1/N,N))
  }
}

We decided to recode the EW strategy, but we could have used the weights_ew() function instead. Finally, the main loop is only marginally different from the single strategy loop.

for(t in 1:length(t_oos)){
  temp_data <- returns %>% 
    filter(Date < t_oos[t]) %>%                         # Same first step
    na.omit() %>% 
    select(-Date)
  realised_returns <- returns %>%                       # Same third step: returns
    filter(Date ==  t_oos[t]) %>% 
    select(-Date)
  for(j in 1:nb_port){                                  # This is the novelty: we loop on the strategies 
    portf_weights[t,j,] <- weights_multi(temp_data, j)  # The weights' function is indexed by j
    portf_returns[t,j] <- sum(portf_weights[t,j,] * realised_returns)
  }
}
apply(portf_returns, 2, perf_met) %>%                   # Taking perf metrics
  unlist() %>%                                          # Flattening the list
  matrix(nrow = 2, byrow = T) %>%                       # Ordering them
  `colnames<-`(c("avg_ret", "vol", "SR", "VaR")) %>%    # Adding column names
  data.frame()                                          # Converting to dataframe

apply(portf_weights, 2, turnover_simple, t_oos = t_oos) %>% unlist()

[1] 0.2189094 0.0000000

We recall the order of strategies: MSR first (line) and EW second (line).

Again, turnover should (and will) be added to the performance metric function. Since uniform weights are constant, their simplified turnover is zero. In practice, that is not the case because weights evolve according to asset returns. The adjustment would imply only a small turnover.

There is a well-documented substantial difference between the two weigting schemes in terms of asset rotation. The MSR that we compute is clearly not competitive: even before transaction costs, its Sharpe ratio is smaller than that of the EW portfolio (see DeMiguel et al. (2009) for further evidence on the robustness of the EW portfolio).

5.2 Make do without loops: functional programming

Finally, we show how to bypass the loop over dates. While this is not useful on small datasets, it can save time on large databases because loops are notoriously time-consuming. Below, we show how to proceed with the map() function in the simple case with only one strategy.

# We create a function that will compute returns for each date:
port_map <- function(t_oos, returns, weight_func){ # Weighting function as argument!             
  temp_data <- returns %>% filter(Date < t_oos)    # Still expanding window...
  portf_weights <- temp_data %>% 
    select(-Date) %>% 
    na.omit() %>%
    weight_func()
  realised_returns <- returns %>% 
    filter(Date ==  t_oos) %>% 
    select(-Date)
  return(sum(portf_weights * realised_returns))
}
# the map() function does it all!
portf_returns <- t_oos %>%                         # The variable over which we loop
  map(~port_map(.x, 
                returns = returns, 
                weight_func = weights_msr)         # MSR weights
  ) %>%                                            # Is sent to the map() function
  unlist()                                         # The output is flattened
perf_met(portf_returns)                            # Compute the perf metrics

6 Characteristics-based portfolio choice

In this section, we start the chapter of strategies based on firm characteristics (features). As an illustration, we check a well-documented (though still controversial) anomaly: the size effect. We build two portfolios: in the first (resp. second) one, we invest in the firms that have a below (resp. above) median market capitalisation. The first portfolio will be a ‘small’ portfolio and the second one a ‘large’ one. Stocks are equally-weighted inside the portfolios.

First, we prepare the variables and define the weight function.

nb_port <- 2                                                      # Nb of portfolios
portf_weights <- array(0, dim = c(Tt-1, nb_port, length(tick)))   # Where we store weights
portf_returns <- matrix(0, nrow = Tt-1, ncol = nb_port)           # Where we store returns

weights_cap <- function(data,j){    
  # More general than before: we feed all the data, not just returns
  m <- median(data$Mkt_Cap)       # Compute the median market cap
  n <- nrow(data)                 # Compute the number of assets
  if(j == 1){return((data$Mkt_Cap < m)/n*2)} # Small cap
  if(j == 2){return((data$Mkt_Cap > m)/n*2)} # Large cap
}

There will only be Tt-1 dates since we lose one because of the computation of future returns.

Second, we can launch the backtesting loop.

for(t in 2:length(t_oos)){
  temp_data <- data %>% filter(Date == t_oos[t-1])        # We keep the data of the previous date
  for(j in 1:nb_port){                                    # We loop on the strategies 
    portf_weights[t-1,j,] <- weights_cap(temp_data, j)    # The weights' function is indexed by j
    realised_returns <- returns %>% 
      filter(Date ==  t_oos[t]) %>%                       # Shift in the date
      na.omit() %>% 
      select(-Date)
    portf_returns[t-1,j] <- sum(portf_weights[t-1,j,] * realised_returns)
  }
}

Third, we proceed to performance metrics.

apply(portf_returns, 2, perf_met) %>%                   # Taking perf metrics
  unlist() %>%                                          # Flattening the list
  matrix(nrow = 2, byrow = T) %>%                       # Ordering them
  `colnames<-`(c("avg_ret", "vol", "SR", "VaR")) %>%    # Adding column names
  data.frame()                                          # Converting to dataframe

Small firms do indeed generate a higher level of performance! In order to reach this conclusion in a rigourous fashion, we would need to perform the same analysis on at least 1,000 stocks (ideally, more) and on 5 to 10 portfolio sorts (from very small firms to very large ones). Below, we check the weights of the portfolio on one particular date.

small <- portf_weights[2,1,] # t = 2, j = 1 (t_oos[2] = 2010-03-01, small firms)
large <- portf_weights[2,2,] # t = 2, j = 2 (t_oos[2] = 2010-03-01, large firms)
data.frame(small, large, row.names = tick)

Indeed, some stocks have zero weights and others 1/15.

Finally, let’s see how we could have coded those strategies using the tidyverse (and a lot of piping!).

data %>% filter(Date > sep_date) %>%                      # Keep only the out-of-sample backtesting dates
  group_by(Tick) %>%                                      # Group by stock
  mutate(F_Return = dplyr::lead(Return)) %>%              # Compute forward (i.e., realised) return
  na.omit() %>%                                           # Take out NAs
  group_by(Date) %>%                                      # Group by dates
  mutate(Mkt_Cap_Binary = Mkt_Cap < median(Mkt_Cap)) %>%  # Compute median cap for each date
  ungroup() %>%
  group_by(Mkt_Cap_Binary) %>%                            # Group by Mkt_Cap: small vs large
  summarise(avg_return = mean(F_Return))                  # Simple pivot table

It can be useful to see how often stocks switch from one family to another (from below median to above median or vice-versa). Below, we show a plot of Mkt_Cap, conditional on Mkt_Cap being above the current median.

data %>% group_by(Date) %>%                          # Group by date
  mutate(Mkt_Cap_median = median(Mkt_Cap)) %>%       # Compute median cap for each date
  filter(Mkt_Cap > Mkt_Cap_median) %>%               # Keep only the large stocks
  ggplot(aes(x = Date, y = Mkt_Cap, color = Tick)) + # Plot
  geom_line() + ylim(75000,250000) +
  geom_line(aes(x = Date, y = Mkt_Cap_median), color = "black") 

# The black line shows the running median

The straight lines show the discontinuities: one stock being large at some point in time, then small and then large again. The straight lines show the periods when the stock was small. The black line shows the median capitalisation (in the sample). Finally, because we focus in the zone close to the median and impose an upper limit of 250B$, there are some missing points.

7 Sustainable investing

ESG investing has gained a lot of traction (see my review paper: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3715753). In this section, we briefly explore this topic.

First, a glance at the data. Let’s look at the distribution of ESG_rank (from Sustainalytics).

load("data_esg.RData") # Loading the data
data_esg %>%
  mutate(Year = year(Date)) %>%
  ggplot(aes(x = ESG_rank, fill = as.factor(Year))) + 
  geom_histogram() + labs(fill = "Year") + theme_light()

There are few firms with extreme ESG_rank (below 10 or above 90).

Let’s continue with a “simple” piping flow that creates a compact dataset.

data_esg2 <- data_esg %>%
  select(Tick, Date, Close, ESG_rank) %>%            # Keep relevant columns
  filter(Date > "2014-02-01") %>%                    # Data available after 2014-02-01
  group_by(Tick) %>%
  mutate(Return = dplyr::lead(Close)/Close - 1) %>%  # Forward return
  ungroup() %>%
  group_by(Date) %>%                                 # For each date: create ESG groups
  mutate(ESG_group = if_else(ESG_rank > median(ESG_rank), "High_ESG", "Low_ESG")) %>%
  ungroup()
data_esg2

Based on this short dataset, let’s have a quick look at the perf of the 2 groups (high ESG versus low ESG).

data_esg2 %>%
  group_by(Date, ESG_group) %>%
  summarise(Return = mean(Return, na.rm = T)) %>%
  group_by(ESG_group) %>%
  summarise(avg_return = mean(Return, na.rm = T),  # Be careful with na.rm = T
            vol = sd(Return, na.rm = T),
            SR = avg_return/vol)

Two conclusions:
1. High ESG does slightly better.
2. But the difference is not huge (t-stat => exercise!)

From individual annual returns: is there are pattern?

data_esg2 %>%
  mutate(Year = year(Date)) %>%
  group_by(Year, Tick) %>%
  summarise(avg_ESG = mean(ESG_rank, na.rm = T),
            avg_return = mean(Return, na.rm = T)) %>%
  ggplot(aes(x = avg_ESG, y = avg_return)) + geom_point(alpha = 0.3) +
  geom_smooth(method = "lm") + theme_light()

No pattern at all, it seems…

Let’s dig deeper and have a look at screening/selection intensity. We compute the average returns of firms that have ESG_rank above a given threshold. To do that, let’s make a function.

library(scales) # This will be used in the plot below to zoom on the values with special limits
ESG_impact <- function(level, data_esg2){
  data_esg2 %>%
    filter(ESG_rank > level) %>%   # Filter to keep only the firms with ESG > level
    pull(Return) %>%
    mean(na.rm = T)
}

level_esg <- 10*(1:9)                                 # ESG_rank threshold
map_dbl(level_esg, ESG_impact, data_esg2 = data_esg2) # The raw result

[1] 0.01106186 0.01111554 0.01128353 0.01122823 0.01153817 0.01141918 0.01184463 0.01164315 0.01184526

level_esg %>%
  bind_cols(return = map_dbl(level_esg, ESG_impact, data_esg2 = data_esg2)) %>%
  ggplot(aes(x = level_esg, y = return)) + geom_col(alpha = 0.5) + 
  theme_light() + scale_y_continuous(limits=c(0.008,0.012), oob = rescale_none)

Overall, higher ESG selection seems promising indeed (but let’s not jump to general conclusions). Finally, let’s have a look at more homogeneous groups. We need different reference points (level_esg) because of the distribution of ESG_rank.

ESG_quant <- function(level, width, data_esg2){ # A new function
  data_esg2 %>%
    filter(ESG_rank > level - width/2,
           ESG_rank < level + width/2) %>%     # Filter to keep firms with ESG around level
    pull(Return) %>%
    mean(na.rm = T)
}
level_esg <- quantile(data_esg$ESG_rank, seq(0.1,0.9, by = 0.1))
level_esg %>%
  bind_cols(return = map_dbl(level_esg, ESG_quant, width = 10, data_esg2 = data_esg2)) %>%
  ggplot(aes(x = level_esg, y = return)) + geom_col(alpha = 0.5) + 
  theme_light() + scale_y_continuous(limits=c(0.007,0.0125), oob = rescale_none)

New names:
* `` -> ...1

The effect is not monotonic: things are not that simple…

8 Exercises

8.1 Rolling window

Change the main loop so that only 60 points of data are given to the weighting scheme(s). Sixty points amount to 5 years of monthly data.

8.2 Other performance metrics

Using the PerformanceAnalytics package (install it first!), compute the maximum drawdown via: https://www.rdocumentation.org/packages/PerformanceAnalytics/versions/1.5.2/topics/maxDrawdown

Compute the transaction cost-adjusted SR: \(TC-SR=(\bar{r}-0.005*Turn)/\sigma\).

Add both metrics to the perf_met() function.

8.3 Minimum variance

Add the MV portfolio to the set of strategies. The weights depend only on the covariance matrix: \(w=\frac{\Sigma^{-1}1}{1'\Sigma^{-1}1}\).

8.4 map() expertise

Extend the map() syntax to the case with many strategies. BONUS: have a look at pmap()!

8.5 Realistic portfolios

In practice, many saveguards are applied, if only to reduce turnover. One such example is box constraint: the weights in the portfolio must not lie above or below user-specified thresholds.

In minimum variance, it is possible to reduce leverage by shrinking the covariance matrix towards the identity functions: \[\Sigma^*=(1-\alpha) \Sigma^S+\alpha I,\] where \(\Sigma^S\) is the sample covariance matrix and \(I\) the identity matrix. As \(\alpha\) increases to 1, the weights converge to the equally-weighted portfolio.